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SYNOPSIS 


In the present work, buoyancy-driven motion o-f 
deformable bubble through a liquid metal bath during co 
conversion has been studied. The solution involves viscous 
past deformable bubbles . The effect of volume of the bubble, 
its position in the bath T on the shape has been investigated, 
shape of the bubble and the flow around it is dealt in a cou 
manner . During different stages of copper making the reac 
species are transported to the bubble— 1 i quid interface, 
unsteady transport of the species and its distribution around 
bubble has been studied. Since no suitable reaction constant 
available in literature at this stage , it is assumed that 
reaction is mass transfer controlled. 

The analysis involves solution of Navier— Stc 

equations , Continuity equation , Laplace equation and diffus 
equations . This also involves certain relationships concerning 
differential geometry of surface of revolution . 

For handling such moving boundary problems. Fin 
Element Method has been used. During the solution each time 
geometry changes , and hence, the problem calls for an automa 
mesh generation and coordinate fixing . For this pur p 
transf inite interpolation is used. 

Since much experimental information is not available > 
such a system to check the validity of the formulation i 
calculation schemes , prediction of a rain drop shape has b< 
compared with experimental photograph which matched to a very hi 
degree . Some limiting deformations has been considered. This 
because of the formidable dif f i cul t ies encountered in solving t 
Navier-Stokes equations for flow around highly distorted bubbles 



CHAPTER I 

INTRODUCTION 


Bubbles play a key role in many natural phenomena and in 
a host of industrial applications. Reacting gas bubbles are of 
fundamental importance in many chemical and metallurgical 
operations such as absorption, floatation, degassing, steel making 
and copper conversion etc. There is relative motion between the 
bubbles and the surrounding fluid in all these cases. 
Consequently, the phenomenon of interest is the transfer of heat 
or species from one side of the bubble interface to the other. In 
the present study the interaction of bubbles in liquid matte in a 
copper converter is dealt with. The hydrodynamics of bubble motion 
and its bearing upon the rate of oxygen depletion from one bubble 
are highlighted. 

1.1 Copper Conversion 

The copper converting operation is carried out in a 
horizontal basic refractory lined Pier ce-Smith type side blown 
converter. The matte that is charged into the converter consists 
of FeS and Cu^S. Converting operation removes the iron and sulfur 
from the matte and results in the production of a crude blister 
copper. The converting operation consists of two stages namely: 

(a) Slag forming ^Th i s is the stage in which FeS elimination 
occurs. The probable reaction is 

2 FeS + 3 O + SiO = a FeO. SiO. + a SCT (1.1) 

2 ' 2 2 2 

(in matte) air) (flux) (slag) (gas) 

In this stage FeS is oxidized to FeO and Fe^O^ and sulphur is 
removed as SO £ . Silica flux is added for combining with FeO (and 
partly with Fe^O^) for slag formation. Form thermodynamic 
considerations it can be shown that this stage ceases when the 



matte contains less than 1% FeS. Liquid fayalite slag becomes 
saturated with magnetite and the principal product of the slag 
forming stage is the "White Metal" or impure liquid Cu^S. The 
remaining sulfur is oxidized to SO^ during copper making stage. 
At the end of the slag forming stage, the converter contains 
molten slag or fayalite (2FeO. SiO^) with 10—20 V. solid 
magnetite (saturated) and up to 15 7. dissolved and entrained copper 
and molten matte which is principally Cu^S, with less than 1 .554 
FeS C76Bis3. These two phases are immiscible with the slag 
floating on top of the matte. 

(b) Copper making s Metallic copper is formed by the 

following reaction. 

Cu S C1D + O.Cg) =2Cu C13 + SO Cg> (1.2) 

2 2 2 

Initially, when air is blown through Cu^S, sulphur is removed as 
SO,, to give a sulphur— deficient white metal. No metallic copper 
is formed during this stage. The overall reaction is as follows: 

Cu £ S + X0 2 = CUg S 1 _ }< + X S0 2 (1.3) 

Subsequent blowing of air causes the appearance of a second liquid 
phase of blister copper containing 1-2 V. S . The metallic blister 
copper phase is more dense than sulphur deficient Cu 2 S phase and 
sinks to the bottom. Further blowing of air results in additional 
sulphur being removed from the system. Eventually the system 
becomes so sulphur deficient that the sulphide phase disappears 
and only blister copper remains. The final sulphur is removed by 
further blowing and the converting process is carried out till a 
trace of Cu 2 0 appears. 

l«2 Formation and motion of bubbles in submerged blown systems 

Considerable amount of both theoretical and 

experimental literature is available regarding the formation of 
bubbles from submerged nozzles. Indeed, this has become a separate 
Field of study. When liquid does not wet the nozzle, the bubbles 
tend to form at the outer circumference of the nozzle. They grow 



up to a critical size and then get detached from the nozzle, 
critical C79SzeH size is given by the relation! 
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where d is the outside diameter of the nozzle 
n,o 


The turbulence established during the blow helps in 
detaching the bubbles. In metallic systems a higher surface 
tension prevails at the gas-liquid interface which tends to 
prevent the break up of bubble thus large size bubbles are 
formed in such situations. 


After getting detached from the tip of tuyuers, the 
bubbles rise towards the free surface due to buoyancy. During the 
rise the large bubbles often get fragmented into smaller bubbles 
and due to buoyancy they rise up. The velocity with which the 
bubbles rise through the liquid is principally determined by the 
buoyancy force that drives the bubble upwards and the viscous and 
form drags which retard the bubble motion. When these 
counter— acting forces balance each other the bubble rises with a 
steady velocity. The important fact associated with the motion 
are! 

(a) The bubble is not rigid and thus the forces acting on it may 
deform its shape. 

(b) The gas contained in the bubble under goes circulation which 
in turn would affect the drag force. 

(c) The volume of the bubble increases as the bubble rises towards 
the free surface. This, in turn, influences the drag phenomena and 
the bubble velocity. 

The above mentioned factors affect the system 
drastically. The deformation becomes severe in a high density 
system like metallic bath and consequently influences the flow 
field and transport of species and heat to a large extent. The 
flow in general, can be characterized by the dimensionless groups 
of Reynolds number, Weber number, BtvBs number and Horton number. 
The deformation phenomena and their dependence on the 





characteristic numbers will be discussed in more detail in latter 
chapters, as these form major objectives of the present study. 

The drag force affects the motion of the bubble to a large extent. 
The total drag coefficient for solely moving solid spheres given 
by the expression 

= 24/ Re (1.5) 

DSt 

according to Stokes law. Two third of this drag arises from skin 
friction and one third is due to form drag. The corresponding 
terminal velocity is given as: C78Cli3 

2 

U JS = gd A p/ 18 ix (1.6) 

The Hardmard-Rybczynski C78Cli3 theory predicts that the 
terminal velocity of a spherical bubble should be about 50 X 
higher than that of a rigid sphere of same size and density. 
However it is commonly observed that small bubbles and drops tend 
to obey Stokes law rather than the corresponding 
Hardmard-Rybczynski results. Moreover, internal circulation is 
essentially absent. Three different mechanisms have been proposed 
for this phenomenon. 

Bond and Newton C7BCli3 found that small bubbles and drops 
followed Stoke's law while, with increasing diameter, there was a 
rather sharp increase in velocity towards the Hardmard-Rybczynski 
value, they suggested that a circulating particle requires energy 
locally to stretch on inter— facial area element over the leading 
hemisphere, which subsequently shrinks over the rear hemisphere. 
It was hypothesized that this process caused additional tangential 
stresses to retard the particle. 

Boussinesq C78Cli3 proposed that lack of internal 
circulation in bubbles and drops is due to an inter— facial 
mono-layer which acts as viscous membrane. A constitutive equation 
involving the two parameters of surface shear viscosity and 
surface dilational viscosity was proposed for the interface, in 
addition to the surface tension. 



The most reasonable explanation for the absence of 
internal circulation for small bubbles and drops was provided by 
Frumkin and Levich C78CliD. Surface active substances tend to 
accumulate at the interface between two fluids, thereby reducing 
the surface tension. When a drop or bubble moves through a 
continuous medium, it adsorbs surface active materials at the 
interface. The adsorbed particles are swept to the rear leaving 
the frontal region relatively less contaminated. The resulting 
concentration gradient leads to a tangential gradient of surface 
tension which in turn causes a tangential stress, tending to 
retard surface motion. 


1. 3. Dispersed Bubble Systems : 

Isolated bubbles behave in quite a different manner from 
a dispersed system of bubbles. In the case of an array of bubbles, 
each bubble will be affected by the influence of its neighbors. 
The behavior of such systems is quite complex and not well 
understood. Most of the experimental measurements are restricted 
to the water-air system and that too essentially at room 
temperature. Extrapolation of such findings to a high temperature 
metallic system is highly impracticable. 


Dispersed bubble systems are usually formed by the 
injection of a gas stream into a liquid through one or multiple 
orifices. In some cases however bubbles are formed through 
heterogeneous nucleation. The principal parameters used for the 
characterization of the multiple bubble systems are the 
superficial velocity, the mean velocity, the void fraction or gas 
hold up, the surface area per unit volume and the mean bubble 

vol. flow rate of gas 
cross sectional area of vessel 

Vol. of gas bubbles 
Vol . of gas and liquids. 


diameter. 
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BS 


and 


where 
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BS 


superficial velocity. 
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g 


void fraction. 


£ 


Dispersed bubble systems may be divided somewhat 
arbitrarily ass — (i) Bubbling Systems (.£ <0.3), (ii) Froths (.£ =0.4 

g g 

- 0.6), and (iii) Foams {£r g = " 0.98). Among metal processing 

systems with argon stirred ladles, the gaseous deoxidation of 
copper and copper converting system will fall into category (i) i. 
e. the bubbling systems. 

1. 4. Literature Reviews- 

Most of the experimental and theoretical studies of 
bubble formation and chemical reaction at the bubble-liquid 
interface have appeared in chemical engineering and fluid 
mechanics literature. Many of these studies have been concerned 
with aqueous and organic liquids. The simple theoretical model of 
Davidson and co— workers was found to describe bubble formation in 
mercury at room temperatures CSIAshT. However when applied to 
industrial converters the model predicts a bubble frequency S to 3 
times lower than that measured. Since Davidson model can only be 
correctly applied to a stagnant, isothermal and non— reactive 
system, this discrepancy is not surprising. 

Ashman and Co-workers C81AshU studied bubble formation 
in a copper converter. According to them, bubbles that form at 
tuyuers of a copper converter are large and their size is strongly 
dependant on the volume rate of air flow, heat transfer and bath 
circulation velocity. One important aspect of these findings is 
that the bubbles occupy a significant fraction of the liquid space 
between the tuyuers and the bath surface. About 40 X of oxygen is 
consumed during bubble formation at the tuyuers and the rest reacts 
during rise and break-up of the bubble. 

Most of the experimental as well as the theoretical 
results available are for the water-air system and theses have 
been usually represented in terms of the following dimensionless 
groups. 

p d U 

Reynolds Number, R = Cl. 7) 

A 
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In the above equations, the characteristic length d =(6 

1/2 6 
v/rc ) is the volume equivalent diameter of a bubble of volume V 

and U is its terminal rise ; p v and <t are respectively, the 

density viscosity and surface tension of the liquid and g is the 

acceleration due to gravity. Fig-1.1, gives a pictorial summary of 

the shapes of the bubbles in different regimes of these 

non-dimensional groups. 

Till date, the number of theoretical papers dealing with 
ideal flow about bubbles is large and the viscous flow past 
deformable bubbles has seldom been investigated. This is 
connected with the formidable difficulties that are encountered in 
solving the Navier— Stokes equations with deformable free 
boundary. The class of problems with a deformable boundary is 
characterized by the existence of one or more boundaries of the 
flow domain whose shape is dependent upon the viscous and pressure 
forces generated by the fluid motion. In this case the shape of 

the boundary and the form of the velocity and pressure fields in 

the fluid are intimately connected and one must solve for the 
boundary shape as a part of the overall solution of a particular 
problem. The most common problems of this type in fluid mechanics 
occur in the motion of two immiscible fluids which are contiguous 
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Fig. 1.1. Shape Regimes for Bubbles and Drops in 
Unhindered Gravitational Motion through 
Liquids. 




with a common interface e. g. gas— liquid interface. For this 
reason a number of limiting cases have been treated theoretically 
either at low or very high Reynolds numbers in various 
combinations with or without def ormability of the interface? but 
no satisfactory theory exists as yet, except for the cases of very 
small deformations. 

The existing published literature on free boundary 
problems in fluid mechanics is quite extensive in number, but 
limited in scope. Three distinct solution methods can be 

identified. By far the majority of papers are concerned with 
asymptotic solutions or limiting solutions in which the interface 
shape while unknown, deviates only slightly from some predefined 
conf iguration. In the case of bubble motion, a number of authors 
have used the so called "domain perturbat ion" method to solve for 
the first deviation from a spherical bubble in a variety of 
flows C34Tay3, C64Acr3. Levich solved the problem of boundary 
layer flow past Spherical non— def ormable bubbles at very high 
Reynolds number and very low Ueber number C59Lev3. Moore 
acknowledged to some extent the deformation of bubble surface, 
considering the boundary layer around an elliptic gas bubbleE65 
Mooli. Elswawi extended Moore's analysis and accounted for the 
interaction between boundary layer and the shape of bubbleE74 
ElsT. In addition a similar approach has been used to consider 
the first deviation from the limiting form of a slender body with 
an arbitrarily small radius-to-length ratio, which is relevant for 
example, to low viscosity drops in uniaxial extensional flows with 
a sufficiently high strain rate C64 Tay3 and C78 AcrZJ. 

A second method of solution suitable for free boundary 
problems is the boundary integral technique which is restricted to 
limiting cases of either zero or very high Reynolds number, where 
the governing differential equation reduces to Laplace's equation. 
In this method, fundamental solutions of the linear governing 
equations are used to reduce the general n-dimensional problems to 
the solution of (n— 1)— dimensional integral equations. The boundary 
integral method is not restricted to small deformations. Indeed 
solutions have been obtained which exhibit large departures from 
predefined shapes C76You3, C81Mik3, C82Lee3. In the analysis 
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done by Miksis, Vanden-Broeck and Keller, they have considered 
very high Reynolds number flow. The effect of viscosity are 
confined to a thin boundary layer around the bubble. Then, by 
treating the surrounding fluid to be incompressible, the motion can 
be represented as potential flow which automatically satisfies the 
Navier— Stokes equation. One way to complete the formulation is to 
match this potential flow to the flow in the boundary layer. But 
it is difficult to do so because even for a spherical bubble the 
flow in the boundary layer is not known near the rear stagnation 
point on the axis of symmetry at 5 = n. As an alternative they 
have assumed the potential flow right up to the boundary which 
yields an exact solution of the Navier-Stokes equations throughout 
the fluid, with vanishing normal velocity on the surface of the 
bubble. However the restriction to potential flow reduces the 
usefulness of the study and for non— spherical bubble the potential 
flow assumption is not appropriate. 

The third and most important class of free-boundary 
problems, which are not restricted to any limiting cases has been 
studied recently. Ryskin and Leal used finite difference method 
which incorporates a numerically generated orthogonal co-ordinate 
system, which is "boundary-f itted" in the sense that all the 
boundary surfaces of the solution domain (including the free 
boundary whose shape is determined as a part of the solution) 
coincides with a co-ordinate line (or surface ) of the co-ordinate 
system C84Rys3. Thus the problem of interpolation between nodal 
points of the finite difference grid when the latter is not 
coincident with the physical boundaries is avoided altogether. 
Indeed, the existence of the interpolation problem in the first 
place seems to be a consequence of the common analytically 
generated co-ordinate systems, such as cylindrical and spherical 
systems, when the latter do not correspond to the natural 
boundaries of the solution domain. The mapping procedure was 
restricted to two dimensional and axisymmetric flow domains. The 
boundary fitted co-ordinate system is denoted as (£ , 17 , 4> > with <fr 
symmetry. In view of the assumed ax i symmetry, these boundary 
fitted coordinates can be connected with the common cylindrical 
coordinates (x, r, via a pair of mapping functions x( ? , 7 ) )and 
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<y( Z , rj) which satisfy the Laplace equations in eovarient form! 


a 

r y 

ax 1 

*JL\ 

r 1 

ax 1 

az 

t f 

sz J 

an 1 

[ / 

an j 

a 

f , 

a& 

. a 1 

r 1 

a# 1 

az 

l y 

az J 

an 

17 

an j 


( 1 - 12 ) 

(1.13) 


where /( ?» 17 ) is called the distortion function 

With a similar approach Christov and Volkov C85Chrl , considered 
viscous flow in intermediate Reynolds number and for relatively 
low Weber number. The problem of the moving boundary is tackled 
by means of an appropriate coordinate transformation. They 
considered a form of coupled boundary value problem for estimating 
the stream function, vorticity function and the shape of the 
bubble. The unknown shape which is defined as r = R {& ) poses 

considerable difficulties for which they introduced the 
transformation of independent variables by the expressions r = n 
R ( & '},&=& . Thus the problem was transformed to a simple 
computational domain for solution purposes. 

Though recently the steady rise of bubbles in viscous 
liquids has been intensively studied both theoretically and 
experimentally, the knowledge gathered about the shape of bubbles 
and patterns of the flow is still not complete. Experimental 
studies devoted to the motion of deformable drops or bubbles have 
been chiefly concerned with the estimation of the terminal 
velocity of rise . Only recently , some quantitative information 
concerning the shape of bubbles and the streamlines of flow 
around them, has been experimentally collected by Bhaga and Weber 
ESOBhaH . They have presented experimental data in the form 
correlations for bubble rise velocity and shape. The basic 
techniques consist of photographing a rising bubble with a cine 
camera moving upwards at the same speed as the bubble. The rate of 
•scent of the camera could be varied. The camera elevator speed 
was set equal to the subject bubble speed from earlier calibration 
runs. Using this technique the bubble remained in the center of 
the frames in the cine film throughout its rise. The velocity 


field around the rising bubble was visualized by hydrogen bubbles 
introduced electrolyti cal ly into the liquid. The liquid used was 
aqueous sugar solution with differing concentrations. 

But much information is not yet reported about the shape 
of bubble and nature of flow in a liquid metal bath. There are 
evidences that bubbles in liquid metal show the behavior expected 
from more concentrated liquids. Because of large surface tension 
forces for liquid metals, the Morton number tends to be low. 
However metallic systems are prone to contaminations by 
surface active impurities which in turn, alter the bubble shapes 
drastically. The picture of two dimensional bubble in liquid 
mercury has been presented in fig-1.2 C78Cli3. For experimental 
convenience the bubble studies have usually dealt with rather 
large bubbles and , there is only limited data available for 
spherical or slightly deformed ellipsoidal bubbles in liquid 
metals. These data tend to suffer from abnormally high scatter 
due to experimental difficulties associated with opaque media and 
high temperature. Moreover, the cold model experiments are not 
solely concerned with steady rise, since during motion the bubble 
passes to regimes of lower pressure and hence the bubble volume 
increases. This is quite significant in molten metal baths of 
high density. The expansion would of course be greater if the melt 
is at a reduced pressure, as in the case of a vacuum degassing 
system. 


Theoretical studies regarding spherical cap bubbles are 
still many. The well known analysis of the “spherical cap” bubble 
due to Davis and Tayler C50Dav3 relates the speed of rise to the 
radius of curvature of the bubble at front stagnation point C66 
Coin, but the overall spherical cap shape is assumed a priori, 
rather than being determined as apart of full solution. According 
to Collins a gas bubble in liquid is termed large if the Reynolds 
number and Weber number associated with its motion are 
sufficiently high so that viscous and surface tension forces are 
insignificant compared with the inertial and gravity forces in 
determining its shape and velocity. The bubble then adopts the 
characteristics of spherical cap form discussed by numerous 
authors. In his theoretical study, Collins used a formerly 


or 



Fig. 1.2. Two-dimensional Nitrogen Bubble 
in Liquid Mercury. 




developed equation which relates the bubble velocity U and the 
apparent curvature a of the spherical cap, in the form : 

U = 2/3 (g a ) 1/2 Cl. 14) 

Richardson and Bradshaw as quoted in C67Rip3, developed 
an equation for describing the mass transfer to a gas bubble from 
molten metals at reduced pressures. Here it was assumed that the 
pressure inside the bubble was identical to the pressure inside 
the adjacent liquid. This is unlikely to be valid in the immediate 
vicinity of the free surface held at a low pressure. It was found 
that significant distortion of the bubble occurred on approaching 
the free surface; further more the observed expansion was much 
less than what one could expect from hydrostatic considerations. 
In their physical model, the authors assumed that bubble rises 
with a constant velocity and it maintains sphericity throughout 
the rise even after expansion. 

Chao C69ChaIl studied the transient heat and mass 
transfer to a translating fluid sphere at high Reynolds number. 
Major assumptions used in that analysis were irrotational flow 
field, constant properties, and thin boundary layers. From a 
similarity analysis, the average, outside Nusselt number at 
steady state was shown to be 

i /? i /o 

Nu , = (2/ rc > Pe (1.15) 

1, ss 

while studying heat transfer between a vapor bubble in motion 
and the liquid, Ruckenstein also obtained identical results as 
Chao. Chao found that for sufficiently large Reynolds and Peclet 
numbers the transient response behavior of the thermal/species 
boundary layers of a circulating fluid sphere depends on a single 
parameter Uc/R and is independent of fluid properties. This is in 
contrast to that for the laminar boundary layers at the front 
stagnation of a solid sphere. He concludes that the growth of the 
boundary layer inside the fluid sphere is dominated by the 
diffusion process and is thus governed by the conduction 
transients. For mass transfer also an exactly similar approach was 



used with heat flux substituted by mass flux and Nusselt number by 
Scherood number. 

Levich E62Lev3 examined the steady mass extraction rate 
from a spherical droplet of radius R falling at a constant 
velocity U in another liquid in the creeping flow regime. The 
velocity fields both inside and outside were those of Hardmard- 
Rybczynski C78CliD. The Peclet numbers for mass transfer were 
assumed to be large and the transfer resistances on both sides of 
the droplet surface were comparable in magnitude. Levich had also 
analyzed the transient response behavior of the concentration 
boundary layers of the translating droplet and stated that the 
time required to attain steady state is C1+ fj 3 R/2 U where yt is 
the ratio of the inside and outside viscosities. 

I h 

” iJ i. S. SCOPE OF IHE PRESENT STUDY s 

1 ’ In the present work, buoyancy-driven motion of a 

deformable bubble through a liquid metal bath during copper 
conversion has been studied. In addition to being an inherently 
ms interesting physical problem, the motion of deformable bubbles 
iM also represents a good example of the important class of free 
>li boundary problems in fluid mechanics, from which one may develop 
» a better understanding of both the solution methods and the 
factors that control the boundary shape. The shape is dependent 
upon the viscous and pressure forces generated by the fluid 
motion. In this case, the shape of the boundary and the form of 
velocity and pressure fields in the fluid are intimately 
;* connected, and the boundary shape has been solved as a part of the 
overall solution of the problem. The flow field has been studied 
i both inside and outside. The Finite element and finite difference 
methods have been used in combination for the analysis of flow and 
l heat transfer in both phases. 

While earlier studies used domain transformation for 
, solving the shape, here simple finite difference technique has 
been used to solve for the bubble shape. The changes in the 
bubble volume and the translational velocities during the rise of 
the bubble are taken into account. The volume change is calculated 
from the prevailing hydrostatic pressure of the metallic bath at 


any depth, assuming the entire bath to be isothermal. The change 
in the bubble velocity is determined by incorporating an 
instantaneous balance between the drag and buoyancy force acting 
on the bubble. 

H The finite element grid generation for this problem has 

been handled by the transfinite interpolation scheme. Though the 
study has been done for a single bubble it can be extended for 
1 multi- bubble system also. In addition to the analysis of shape 
¥ and flow field, the mass transfer of reacting species both inside 
and outside has also been studied. Though the usual bubble shapes 
in metallic baths are of spherical cap geometry, it has not been 
practicable to simulate these conditions numerically when the 
when the shape distortion is severe. However, the simulation has 
been successfully carried out for considerably distorted 
1 ellipsoidal bubbles with an equivalent diameter less than or equal 
to 5 mm. 



CHAPTER-II 


MODEL FORMULATION 


It has been experimentally proved that in different 
regimes of Reynolds number and Etttvtts number, the bubble takes 
different shapes. But the basic question is why should the bubble 
have any resemblance to a spherical shape at all ? This is 
because of the dominant influence of surface tension forces, for 
small size bubbles. Surface tension always tend to reduce the 
surface area of free mass of fluid to the smallest value 
possible. The smallest possible surface area for a given volume is 
that of a sphere, and an isolated bubble or drop which is not 
distorted by external forces is pulled by uniform surface tension 
into a spherical shape. In terms of thermodynamics, the interface 
adjusts itself to the spherical shape in order to minimize its 
surface free energy. 

This question also can be considered in terms of 
pressure, i. e. the pressure field prevailing in the surrounding 
fluid. The smaller the radius, the greater is the pressure 

difference between the two sides of the convex surface. For 

example, a cloud droplet of one micron radius has an internal 
pressure of more than two atmospheres. If the isolated bubble 
should momentarily assume some shape other than that of a sphere, 
its surface would have different radii of curvature at different 
points and the internal pressure just below the surface would 

momentarily be dissimilar at these various points. The consequent 
pressure gradients within the bubble would tend to force air from 
the regions of sharp surface curvature. This is equivalent to 
saying that surface tension, through its control of internal 

pressure, reshapes the bubble into a sphere whenever it happens 
to become slightly deformed. When the bubble is finally brought 
into the shape of a perfect sphere, the uniform surface curvature 
makes the pressure difference uniform at all points of its 
surface, and the internal pressure within the bubble is also 
uniform, provided the pressure field around the bubble remains so. 



The disturbing effects on bubble shapes that appear, as one 
considers larger and larger bubbles , seem to be almost entirely 
due to dynamic and gravitational forces. 

It is well known that when an object moves in a fluid 
the pressure just in front of the object becomes higher than 
average and the pressure in the rear region lower than the 
average. This means of course, that the internal pressure within a 
bubble or drop must change accordingly, the bubble develops an 
excess of pressure near its bottom and a deficiency of pressure 
all around its waist. And if one momentarily assumes that the 
fluid flow is that of a perfect fluid, there will be an excess of 
pressure near the top of the bubble as well as near the bottom. 

The mechanism by which the shape of the bubble gets 
suitably adjusted is as follows. The gradients of internal 
pressure drive air from near the base and top out into the regions 
around the waist, thereby tending to flatten the bubble with an 
increase in the horizontal diameter. Even more intriguing is the 
fact that the resultant modification of the bubble surface 
curvature is just of the right kind required to help the bubble 
restore a uniform internal pressure and achieve an equilibrium. 
The sharpening of the curvature around the waist adjusts the 
surface tension effects to make up for the deficiency of external 
pressure there, while the flattening of curvature near the base 
and top tend to cancel the effects of the excessive internal 
pressure in those regions. Together, the joint action of surface 
tension and dynamic forces deforms the bubble continuously until 
it reaches a stable internal pressure distribution. 

Now it may seem that a large bubble has truly brought 
itself into complete mechanical equilibrium when its internal 
pressure is uniform. But it must till meet an important demand, 
that of the hydrostatic head, that is the bubble must develop a 
vertical pressure gradient in such a way as to satisfy the 
familiar hydrostatic equation relating liquid density, liquid 
depth, and acceleration of gravity. This is quite prominent in a 
high density metallic bath as considered in the present study. 



2.1 Assumptions of the model 5 — 

The following simplifying assumptions are made in the development 
of the model. 

(a) One isolated bubble surrounded by the matte is the system 
considered. 

<b) The bubble motion is nearly steady. A coordinate frame (non- 
accelerating) attached to the bubble center is considered. In this 
coordinate frame the bubble is stationary and the fluid is 
flowing past the bubble in the direction opposite to bubble 
motion . 

(c) The flow system possesses axial symmetry. 

(d) The mass diffusivity, density, viscosity and surface tension 
of the gas and the matte are independent of temperature and 
concentration. 

(e) The reaction at the interface is mass transfer controlled. 

(f) At bubble surface, the surface shear stress is zero, because 
of the large deference in viscosities of the two fluids. 

(g) The bath is isothermal , because of high the conductivity of 
the matte. 

(h) The bubble surface is free from impurities. 

(i) There exists a species boundary layer in the matte phase near 
the bubble interface, because of the large peclet number for the 
liquid. 


An isolated bubble has been considered here. But 
in a practical situation the system consists of many bubbles. 
Hence the flow field of one bubble may influence the other. Though 
the flow domain has been considered to be large, it may be 
compressed by suitable correction for the velocities at the 



established boundary of the flow domain according to the 
celebrated cell model approach. ( for details see appendix— 3.) 

In order to calculate the shape we must consider the 
bubble surface as a moving boundary. In this case the shape of 
the boundary and the velocity and pressure field distributions 
are intimately connected and one must solve the shape as a part 
of the solution. Since the bubble expands at it rises through a 
column of liquid, its motion and shape are not steady. To 
simplify the analysis we have assumed a quasi-steady rise of the 
bubble so that the buoyancy and the drag acting on the bubble 
nearly match each other. From this balance , the bubble velocity 
can be calculated as a function of instantaneous size. For the 
ease of computation, a coordinate frame which moves with the 
bubble is considered. Once the flow is established and shape is 
determined at a height then the distance traveled in that time 
step is predicted and the flow and other calculations are done in 
a quasi— static fashion till the Og in the bubble is exhausted 
completely. 


The mass diffusivity, density, viscosity and surface 
tension change with temperature and the composition of the medium. 
Since much information is not available about the mass 
diffusivity, surface tension and the viscosity data of the system, 
some of the data have been considered in the near order of 

magnitude only . For the effect of change in temperature and 

composition see appendix-2. 

In the slag making stage, FeS is oxidized preferentially 
over CUgS according to the reactions: 

FeS + 1.5 0 2 > FeO + SOg 

6 FeS + 10 0 > 2 Fe„Q„ + 6 S0„ 

c- 3 4 2 

The relative proportion of oxygen consumed and SOg 
produced within the air bubbles are determined by the 

stoichiometry of these reactions and the extent to which each 

occur 

If the composition of the converter slag is known, an 



overall reaction can be obtained from a combination of the above 
two reactions . For a slag containing 23 7 . Fe^O^and 30X FeO , simple 
mass balance gives the overall slag making reaction as s- 

FeS + 1.57 0 2 > FeO + SOg 

since the molar ratio of oxygen to iron is 1.14 in the slag. At 
high temperatures of conversion, the intrinsic rate of this 
reaction will be very fast and hence mass transport is likely to 
be the rate controlling step. Its equilibrium constant is 
extremely large. Therefore the reaction rate must be controlled by 
transport of one or more species, i. e. FeS in the liquid phase 
or Ogin the gas phase. Same arguments hold good for the copper 
making stage also. 

It is easily shown that the transport of FeS to the 
bubble matte interface will be rate controlling step if 


C 


b 

FeS 


< 


1/1.57 



FeS 



Such a condition may be prevailing at the end of the slag 
process. Similarly, for initial stages Ogtransfer within 
bubble will govern the reaction kinetics if 
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b 

FeS 


> 1/1.57 



FeS 



making 

the 


where K = mass transfer co-efficient 
C* 1 = bulk concentration 


Both the transport step will be rate determining under the 
condition that 

Cq (in between) 

The stoichiometry of the blister making reaction indicates that 
certain volume of oxygen removed from the bubble due to reaction 
is replaced by the same volume of SOg. Hence the consumption of 
oxygen has no effect on bubble size during blister making. 

The actual reaction nature has not been reported in 


C FeS = 1/1 - S7 


K, 


2/K 


FeS 
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detail. Since the reaction constant is not available in literature 

the reaction at the interface is assumed to be mass transfer 

controlled. If a suitable reaction constant is made available the 

model will not change drastically and the only changes will be 

incorporation of a source or sink in the mass transport equation. 

For outside it is assumed that there exists a mass boundary layer 

7 

since the mass peclet number is quite high (in order of 10 ) . 
Beyond the boundary layer uniform concentration prevails in the 
matte. 

Traces of surface active contaminants may have profound 
effects on the behavior of the bubble, even if the amount of 

impurity is so small that there is no measurable change in the 
bulk fluid properties. The contaminants cause drastic effects 
such as: 

(i) Retardation of internal internal circulation there by 
rendering the interface rigid. Also the drag is significantly 
increased and the transport rates are decreased. ( ii ) Surfactants 
have greatest influence on terminal velocity near the point of 
transition from rectilinear to oscillating motion. This is 
presumably because internal circulation can drastically alter the 
wake structure of the fluid particle tending to delayed boundary 
layer separation, smaller wakes and delayed vortex shedding. 

(iii) Surfactants play a particularly important role in high 
surface tension systems since surface tension reduction with 
contamination is largest for these systems. 

2. 2. Governing Equations and Boundary Conditions:— 

Consider the quasi-static rise of a bubble in a large 
bath of liquid. Axisymmatry about the direction of motion of the 
bubble has been assumed in this study and cyl indr ical coordinates 
are employed for convenience. 

Since the pressure is a very important parameter for the 
deformation of the bubble the flow solution is obtained directly 
by u, p and v formulation of the Navier— Stokes equations and 
the continuity equation . The governing equations can now be 
written in vectorial form as: 
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and V stands for u and v velocities in Z & r 
directions respectively. 

Though the flow achieves nearly a steady state very 
soon , the species transport is only unsteady in nature . The 
unsteady state mass transport equation in the absence of any 
source or sink term is given by 


3 X 
a "t 
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C V . V X A > 


2 

n> 7 X . (2.3) 
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where X^ stands for the mole fraction of any reacting species i. 
e . Cu^S , FeS or 0^ & H> stands for the dif f usivity of the 
phase considered. For the liquid phase , the dif f usivity is very 
small resulting in a large value of the mass Peclet number . Due 
to this the concentration varies in a thin layer adjacent to the 
bubble surface . In this region , the gradient of concentration 
normal to the bubble surface is much larger than the gradient 
parallel to the surface. The thin concentrat ion boundary layer 
approximation consists of neglecting diffusion parallel to the 
surface in comparison to the diffusion normal to the surface . 
Defining y = r - R, equation 2.3 can be simplified as: 


- a- - t - ■ A + C V . V X.) = ED (2.4) 

at A a 2 

a y 

Formally 7 the above approximation requires Pe — >oo, 

which for most practical situations means Sc kof or any finite 

Reynolds number . Surprisingly this approximation is often 
reasonable E78C1 iT down to Sc of the order of unity. Use of the 
thin concentration boundary layer approximation is sometimes 
called the asymptotic solution for Sc >o&, since it is not 



required that Re be large or that the momentum boundary layer 
approximation be made. 

The representation of the boundary layer adjacent to 
the bubble and the appropriate coordinate system are sketched in 
figure 2.1 . The x— coordinate is parallel to the surface, (x=0 at 
the front stagnation point), while the y— coordinate is normal To 
the surface. The distance from the axis of symmetry to the surface 
is R. It can easily be shown that even for a curved surface , the 
boundary layer equation takes the same form as that of a flat 
plate in terms of the curvilinear coordinates x and y. The final 
form of the mass transfer boundary layer equation is given by 
equation (2.4) . 

These equations are non-dimensionalized using the following 
parameters . 



where 

Re = Reynolds number 
Pe = Mass peclet number 
Fr = Froud number 
t = Characteristic time 
R = radius of the bubble 
w = kinematic viscosity 
U» = uniform velocity 
E6 = Efitttvs number 

R = radius of equivalent spherical bubble 

For the sake of comparison with experimental work, however 
the diameter of the equivalent spherical bubble ( with the same 
volume) is selected for the reference length. The radius of 
equivalent bubble is defined as 




Fig. 2d Coordinates for the thin concentration boundary layer approximation 
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(2.5) 


R = ( 3/4n V ) 

where V = volume of the bubble given by the 
expression 


The non dimensional form of the equations are as follows: 
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The following boundary conditions have been imposed to 
solve the governing equations (2.7 to 2.10) At infinity the 
liquid moves with a uniform velocity unity in the Z- direction. 
Numerical implementation of this condition requires that the 
region be cut at a certain r = r 7 this being of crucial 

<ss 

importance for the efficiency of the cal culation . This distance 
may be compressed by improving the condition in a manner that 
allows a good accuracy even if the solution domain is not very 
large ( see appendix- 3). On the other hand the computational 
efficiency can be improved by generating the grid exponential ly r 
through the transformation r = e" . 

In order to solve the Navier-Stokes equations in 
multi-phase problems » boundary conditions are required between the 
velocities on either side of the interface between the phases . 
The existence of an interface assures that the normal velocity in 
each phase is equal at the interface i . e . 



V = ( V ) C everywhere on interface ) (2.11) 
n n p 

where the subscript n refers to the direction normal to the 
interface. For bubbles of constant shape and size the normal 
velocity is zero with respect to the origin fixed to the bubble. 
The condition of tangential velocity continuity at the interface 
due to no slip condition takes the form 

V t = ( V t ) (everywhere on interface) (2.12) 

X Vr P 

For fluid particles, two additional boundary conditions 
are provided by Newton's third law which requires that normal and 
shearing stress be balanced at the interface separating the two 
fluids. 

The interface between two fluids is in reality a thin 
layer, typically a few molecular dimension thick. The thickness is 
not well defined since physical properties vary continuously from 
the values of one bulk phase to that of the other. In practice 
however, the interface is generally treated as if it were 
infinitesimally thin i. e as if there were a sharp discontinuity 
between two bulk phases. Of special importance is the property of 
inter— facial tension which is best viewed as the surface free 
energy per unit area at constant temperature. 

A complete treatment of inter-facial boundary conditions 
in tensor notation is as follows. If surface viscosities are 
ignored the stress condition reduces to 
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where k is curvature of the surface and t are the deviotoric 

nn 

normal stresses. Under static conditions the above equation 
reduces to Laplace equation. The tangential stress balance 
leads to 
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_ . and t refer to the shearing stresses and 7 
n+ n 3 s 


where t 


is the 



surface gradient. The gas inside the bubble is in— viscid Cp^<< 
and the surface exhibits no intrinsic shear or dialational 
viscosities. Also, usually there is no variation of surface 
tension on the bubble surface if the surface is not very much 
contaminated by surfactants. As a result, the tangential 
components of the stress vector have to be nearly equal to zero 
at the bubble surface. 

t =0 
n+ 

On the same ground equation-2.13 reduces to 

& V 

p -p + 2 m — s — — = & k (2.15) 

r p an 

As already said earlier the air from a high pressure 
region is driven to a region of low pressure in order to maintain 
an equilibrium or uniform thermodynamic pressure inside the 
bubble. Therefore, the internal pressure variation may be ignored 
and the resulting form of Laplace equation becomes 


a v 

-p + E)J . — - — = a k (2.16) 

Cf n 


From differential geometry of surface revolutions , the 
curvature of an arbitrary surface can be defined in terms of 
radius and e in spherical coordinate. (For detail see appendix— 1) 
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( 2 . 17 ) 


Now using this relation, the Laplace equation becomes 
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Non dimensional isat ion of this equation results in 
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This is the fundamental equation which determines the 
shape of the bubble the condition of tangential stress equal to 
zero condition gives the following relation. The tangential stress 
on any arbitrary surface can be defined as 

(o- -t y y 

r =— * ^ sin 2 «+ t cos 2 a (2-20) 

s 2 xy 

(see figure 2.2) 

where & = -p + 2 p Su/Sx and & = -p + 2 u dv/dy 

xx yy r r- j 

r = (j( du/dy + dv/dx ) 

xy 

t Z (du/dx — Sv/dyY) _ . . . _ 

2 — sin 2c* + p ( &u/Sy + 3v/SxJ cos 2c* 

= O 


or 

( Ssf/Sy sin 2c* + dv/dx cos 2a) + ( du./dy cos 2cx - du/dx sin 2*) = O 

( 2 . 21 ) 

where a is the angle between the outward normals at any surface 

The velocity with which the bubble rises up is also a 
function of the drag co— efficient and the buoyancy force on the 
bubble. The drag contribution comes from two sources, namely, the 
viscous drag and the form drag. As the geometry is changing ,this 
force also changes causing a change in the velocity of rise. This 
can be expressed as 

C D = 8/3 (1 - p g/p 1 > gR/ U 2 (2-22) 

The drag force exerted on the bubble surface is of the form 
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where cos at and cos ft are the direction cosines of the normal to 
the surface ds. 

Non dimensional ised form of the expression is as follows 
Fjj/d/2 p U 2 )n R 2 = C D = 

2 £ ^ (— p*+ 2/Re £u?dx?coslx +1/Re (<9u?dy*+dv?dxTcoslyJy*dl 


(2.26) 


where p* = p /pU 


Computing the change in the velocity of bubbl'e rise can be 
predicted. Both and R of the bubble will change during the 
upward rise, and so also the Reynolds number. 


The fact that mass transport is the rate controlling step gives 
rise to the condition that 

X A =° 

at the inte-rface for all time steps. In addition, the mass transfer 
boundary conditions ares 

X A = 1.0 at all other points inside the 

bubble at t = 0 . i. e X ft *= 0 at Y =0 and X* A = 1 as 

y >eo 
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The governing equations and boundary conditions above 
have been solved in a coupled manner using the Finite element 
method. The details of the solution procedure are described in the 
next chapter. 



CHAPTER III 

SOLUTION PROCEDURE 


Analytical solution of the governing equations presented 
in chapter-II is not possible. Hence certain numerical schemes 
have been used. Before solving a physical problem on a digital 
computer, one has to decide which scheme is to be used for 
computation. Variety of numerical techniques are available such as 
Finite difference. Finite element and Finite volume method for 
solving these equations. In the present problem the geometry of 
the physical domain is not regular. In addition to that the 
parameter of importance for shape calculation is pressure and 
velocity. Therefore Finite element method is considered to be the 
most suitable technique for obtaining velocity profile and 
pressure directly for this irregular moving boundary problem. 
While the flow calculation is done using finite element method, 
finite difference method is used for the shape calculation from 
Laplace equation and for the solution of the concentration 
boundary layer in the liquid phase. Therefore the solution 
procedure involves both FEM and FDM in combination. 

3. 1. FEM Analysis^ 

In FEM, the governing differential equations are 
converted into equivalent set of integral equations which are 
minimized over the whole solution domain. The solution domain is 
divided into several small elements of known shape and solution 
variable is approximated within each of these elements by suitable 
approximation of functions. The integral equations are then 
evaluated within each element using the approximating functions 
and all such elemental integrals are assembled to obtain the final 
matrix equations to be solved. Two procedures are available for 
converting the governing partial differential equations for 
equivalent integral equations. These are 

(a) Variational Procedure 

(b) Galerkin's Procedure 



In the present work, Galerkin's approach is used. The 
Galerkin's technique is one among the class of techniques 
collectively known as the method of Weighted residuals C77Zie!3. 
This procedure is used on the fact that, for the differential 
equation of the form 

L C 4 > 1 - O (3.1) 

where 4 > is the exact solution to the differential equation, an 

4f> 

approximation 4 > will leave behind a residue. That is 
L C > R (3.2) 

where R is the residue, which may be a function of spatial 
coordinate and time. In Galerkin’s procedure, the residue weighted 
appropriately by suitable weighting functions is minimized over 
the whole solution domain D. In mathematical form, this becomes 


f W.R dv = f H. L(<£) dv = 0 
J D 1 J D 1 


(3.3) 


for each i , where W. are weighting functions. The approximation 
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function may be written as a piece continuous profile for each 

element as m 
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In the above equation, m is the number of nodes per 

(e) 

element and N. are the interpolating functions or shape 

* * 
functions of the elements. The nodal values of 4 > for the element 

C 6 ) 

are denoted by 4 > - • In Galerkin's procedure, the weighting 

1 ( e ) 

function U. are taken to be the same as the shape functions N. 

l l 

within each element. 


3. 2. Types of formulations* 

There are four types of formulations which 
are in use for obtaining finite element solutions of any flow 
problem. Those ares 

(1) Stream function formulation. 

(2) Stream function and Vorticity formulation. 

(3) Penalty formulation. 

(4) Velocity Pressure formulation. 


The stream function approach for solving the flow 



problem was first proposed by Os 1 onC750s 13 . By introducing the 
stream function , the cont inui ty equation can be satisfied 
exactly and the momentum equation can be combined into a single 
function in terms of y>. The main disadvantage of this formulation 
is that the governing equations and the interpolating functions of 
higher order are required. Also in order to get velocity and 
pressure, further processing is required. 

In stream f unction— vorti city formulation , in stead of 
dealing with a single governing equation of fourth order, the 
stream function (y) and vorticity(w) are chosen as the unknowns and 
two governing equations of second order in y and <•> result 
C73Tay3 . In this formulation, the boundary conditions may consist 
of specified values of up and or specified values of their first 
derivatives. The principal difficulty in using the stream 
f unction-vorticity formulation is that , in general ,vorticity is 
unknown as priori along the boundaries. 

In penalty function C78Hen3 approach, pressure variable 
is eliminated and only the velocity components are solved. This 
approach is sometimes useful because retaining pressure as a 
variable often leads to numerical difficulties. The penalty 
formulation involves treating continuity equation as a constraint 
among the velocity components and the continuity constraint is 
introduced into the momentum equation in a Lagrange multiplier X. 
It is possible to choose X to be very large and eliminate the 
pressure in terms of X. The boundary conditions to be specified 
for this formulation are either specified values of velocity 
components or specified values of their first derivatives. 

In velocity pressure formulation, the pressure and 
velocity components are considered as the nodal variables. The 
advantages of this formulation are : 

( 1 )The formulation is readily extended to three dimension. 

(E)Dnly (linear element) continuity is required for 

element interpolation functions. 

(3>Pressure, velocity, velocity gradient and shear boundary 
conditions can be directly incorporated into the matrix equations. 

(4) Free surface problems are tractable. 

(5) The formulation takes less computational time than the 
stream function vorticity f ormulation. C73TayU 


04 


In the present study, this approach is used because the 
primitive variables are obtained directly. 

3. 3 DERIVATION OF FINITE ELEMENT EQUATION £ 

The discretized domain is shown infig— 3.1a & 3.1b. The 
elements used through out the work is the isoparametric 
quadrilaterals containing eight nodes, one at each corner and one 
at the midpoint of each side ( see fig-3. 2 ). In general terms 

parabolic element is one which can accommodate a parabolic 
variation in the variables along any side of and across the 
element. The term isoparametric implies that both elemental 
geometry and variation in the variables across the elements are 
described using the same type of polynomial or shape function. 
Since the element sides can follow parabolic shapes then quite 
complicated physical geometries can be mapped with few such 
elements . 


A 

5 6 7 

CFig - 3.23 Eight, noded isoparametric element. 

Under most circumstances the utilization of normalized 
coordinates standardizes the integration process. For example 
consider a general form of equation in terms of global coordinates 

J J F’ (x,y)dxdy (3.5) 

a 

or introducing the limits of integration 
x b y b 

•£ 1 F T (x,y )dx dy (3.6) 

a y a 

for one element it can be written in terms of normalized 
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FIG. 3.1b. INSIDE GRID STRUCTURE 
N NODE NUMBER 

(Q) AUTOMATICALLY GENERATED ELEMENT NUMBER 
<fij FIXED SET OF ELEMENT NUMBER 





FIG. 3.1a. EXPONENTIALLY DISCRETIZED GRID 
FOR OUTSIDE. 
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curvilinear coordinates 
+1 +1 

f f FUf,Y>)d?dy> 

-1 -1 

which can be integrated quite readily using simple integration 

techniques. A numerical integration procedure is adopted where the 

sampling points are termed as gauss points. In particular the the 

Gauss— Legendre quadrature is used leading to a high accuracy. For 

example +1 +1 

I = f f F(£ )d£ dr? <3. 8 ) 

- 1-1 


after 


implementing Gaussian 
ngaus ngaus 



integration 
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ngaus = number of gauss points 
a^ & a^ are weighting factors 


and 


are coordinates of the i lh and j* h integration points. 


The shape functions 

Corner nodes: N.= \ ( 1+£ .? ) ( 1+n . 17 ) (? .£ +>? .Y)~1) 
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Mid side nodes: hL = ^ ( IHf 2 ) ( 1 + 17 ^ 17 ) »? = 0 
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Clearly if a local coordinate system is to be used then the 
first order variation with respect to the global coordinates must 
also be expressed in terms of local coordinates as: 
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J=Jacobian evaluated explicitly since the local variation in 
x and y can be defined. 
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dx dy = det J df dr) 

In the present work,the primitive 
variables u, p and v are evaluated from the solution of coupled 
Navier— Stokes and continuity equation iteratively. The procedure 
is then to assume a starting value or initial guess for u, v and p 
then solve for these variables till the solution is converged. 

Following the accepted practice of depicting the 
variation of pressure by shape functions of one order lower than 
those for defining the velocity distribution CSITayH 
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where n = 8 and m = 4 

This requi remen t is effected by utilizing the same element 
geometry where all eight nodes are associated with velocities and 
only corner nodes with pressure. Spatial coordinates within the 
element are defined in terms of all eight nodes . This element is 
called super par ametri c when fewer nodes are used to define a 
variable other than the spatial coordinates . 

Employing the Galerkin's weighted residual approach , eq.2.1 and 
2.2 result in 
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The higher order derivatives require higher order shape 
functions or polynomial s . They can easily be reduced into lower 
order by weak formulation. 
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The resulting equation becomes 
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(3.24) 


If cos ot > sin ot , then (3.24) is u— equation and (3.23) is 
v— equation and vice versa. Implementation of this condition gives 
diagonal dominance to the respective equation. 

The shape deviates from sphericity. Hence the direction 
cosines are now functions of the surface revolution. They are to 
be evaluated on the surface itself. These outward normals are 
calculated in turn and different values are associated with 
different sides as follows. 


y— component of outward normals 
SIDE 
1 
2 

3 

4 

x— component of outward normals 
SIDE 
1 
2 
■O 

4 

The elemental side length are 
locations. These are calculated as 
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the element side making up the boundary of the domain. 


In the present study, the boundary condition have been 
applied directly using eq 2.11 and 2.21 . 
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the expanded form of eq.3.15 to 3. 18 are as follows: 
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The y-component of the momentum equation is 


n e r a 8 3 N. S S 3 N. 

1 [ J fi N i 1 E \ v k E j-?+ E N k Uk E >Vj « 


4 3 U 


+ f N. E M 3 P ■ dV - f N . 


ft 1 


Cl F r 


P 3 N.a 3 N. 3 N . a tfhl. 

+sr L J Q ,r ^ ^ v . jdW * 1 “j^r 0 <3 - E1) 

+ tot 1 

Similarly the x-momentum equation 
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On the boundary 
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8 8 

\ N. cos a u. + V N. sin « v. da = O (3. 24 ) 

4 j j 4 J J 

If cos a. > sin a. , then (3.24) is u-equation and ( 3 . 23 ) is 
v— equation and vice versa. Implementation of this condition gives 
diagonal dominance to the respective equation. 

The shape deviates from sphericity. Hence the direction 
cosines are now functions of the surface revolution. They are to 


be evaluated on the surface itself, 
calculated in turn and different 
different sides as follows . 
y-component of outward normals 
SIDE 
1 
2 

3 

4 

x-component of outward normals 
SIDE 
1 
2 

3 

4 

The elemental side length are 
locations . These are calculated as 
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The final assembled form of the equation is 
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where 1x1 and 1x2 are direction cosines of the x and y global axes 
to the direction of the gravitational field. 

On boundary, the coefficient matrix changes depending upon the 
equation which is being considered as u and v equations. The matrix 
equation is solved iteratively by Frontal method. CFor detail see 
appendix— 4) . 

Once the primitive variables are obtained, the shape of the 
bubble is predicted from the Laplace’s equation. This equation is 
highly nonlinear due to the presence of first and second order 
derivatives of function R(<9). Hence it has to be solved 

iteratively. The construction of the difference approximation was 
really challenging, because of the instable nature of the 

equation. The aim was to obtain as implicit as possible a scheme 
in order to decrease the jump of the solution vector. The 
following quasi linear form was successful to get a converged 
solution of the function r = R(0 ) . One fictitious time step was 

introduced. The simplified form of the equation is as follows. 
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Now implementing the above said implicit scheme the resulting 



orm of the Laplace equation becomes 
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Introducing the difference equation, the resulting equation 
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»akes the form 

P k+1 

R- * +R- „ -2R. -R ) 

i+1 i-1 i i 


( R . „ - R. . > 
i+l 1-1 


+CAO + CR 


R. )* 
1 — 1 


R . — R . 

i+l i—l 


2A& 


2R 


cot & . — R . 

i i 


i+1 


Ap Eo 


2 2 
4A© R. 
l 


(R 


i+1 


8Ad 


R. . ) 2 3/2 
- 1 " 1 - > 3 


(3.27) 


As the function R(£) is to be smooth at the axis of symmetry 
of the spherical coordinates , the boundary conditions for RCd) are 


~-= O at & = O and n 


(3.28) 


which is approximated as 
- 3 Rl - 4R 2 - R 3 = O 


- 3 r n - 4 R N-1 


R N-2 = 0 


so 


that 


the 


equation 


comprise 


tridiagonal 


system which is solved by using the Tridiagonal algorithm. Even 
this solution resulted in a monotonic increment of the solution 
vector for which an additional constraint was imposed on the 
solution for R(S). This is the preservation of the dimensionless 
bubble volume namely 
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from which the characteristic length scale of the problems, the 
equivalent bubble radius is determined as 


.3 V .1/3 

a = (--- ) (3.30) 

4 n 

The monotonic increment of the solution to R(©> was arrested 
by dividing the function R (£) ) by this equivalent radius “a" which is 
determined at each fictitious time step. Simpson's rule was used 
for determination of equivalent volume. Imposition of 


this 


restriction on the solution of R(0) is required because of the 
boundary condition- 

$ R 

Here we have — — — = 0 at & =0 and n . 

This is an over deterministic situation. In such a case the 
solution will involve an arbitrary constant which can not be 
uniquely determined. For example, let us consider an equation of the 
following form 

R ' *= constant 
Now, R *= constant & + C 
But R ’= 0 at © =0 and 
R ' = 0 at & - n 

Implementation of the boundary condition gives the following 
relation 

0 = constant 0 + C 


0 = constant n + C 

Hence C can not be uniquely determined. 

Once the shape and flow field both inside and outside are 
established the mass transport equation becomes linear. 

Applying Galerkin' s approach to eqn 2.9 the resulting 
equation is as 
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Using weak formulation and Green's integration techniques 
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Since there is no gradient boundary conditions, the line 
integral term vanishes, 
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The transient mass transport equation for inside(Gas phase) 
is solved by Implicit Finite element method. The resulting equation is 
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Similarly equation for mass boundary layer is solved where 
the a x term is considered as if the transient term and il 
is marched in x. Here x , is a fictitious time coordinate . 
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where C . = rc/(Pe u Ay ) 
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Drag coefficient is calculated by using the Gaussian 
quadrature formula. Similarly the consumption of oxygen is 



j-stimated by integrating the mass flux along the boundary of the 
mbble . 

The flow solution is done under the assumption that R<0) i s a 
enown function. In the next stage the shape is established by the 
evaluation of the normal stress. However the f l ow field i s strongly 
dependent on the domain of the second phase. Hence shape and fl° tf 
are solved iteratively till both are converged. 

Implementation of such a solution scheme requires automatic 
grid generation. This is done by using ^cansf inite interpolation 
scheme. The details of the interpolation scheme i s described 1° 
Appendix-5. Grid generation using the transfinit e interpolation 
requires that the coordinate of boundary of the doraain be 
specified. For outsider i.e. T for liquid phase, efficiency of the 
computation depends on r = r^ be cut at a considerable distance 
from the bubble but without losing the accuracy. Therefore f° r 
outside, exponential distribution is adopted for elemental siz®' 
The advantages of such a distribution are 

( 1 )The mesh is very fine in the more disturbed region 
, i.e. , adjacent to the bubble surface. 

(2)The number of elements does not increase to a large extent 
so that computational time is less without loss of accuracy. 

However for inside the distribution is linear. Again f° r 
inside the first set of element emerge to be triangular. But two 
triangular elements are combined to form an isoparametric 
rectangular element. This is done just to avoid the complicity that 
would emerge for handling the triangular elements with rectangul ar 
elements adjacent to it. The disadvantages of adoption of such 3 
method are 

(1) The first set of elements that emerges from the center i s 
fixed. 

(2) Very large deformation of the bubble can not be 
considered. 

The grid used is shown in Figure 3.2a and 3.2b.o ne important 
aspect of the grid generation technique is the distribution of 
the first point on both the axes. After the deformation one si<*® 
gets compressed and other side is expanded. Hence for each 
iteration step this has to be adjusted so that the mesh be finest 
adjacent to the bubble surface. 


The numerical calculations have been done using the 
above mentioned formulations for three different bubble 
diameter. The results obtained are presented in the next chapter. 



CHAPTER IV 


RESULTS AND DISCUSSIONS 

Numerical solutions of the model discussed in 
chapter— II were obtained using typical values for the parameters 
of the problem. The physical properties were taken as constants 
for all the computational runs and their values are summarized in 
table 4.1. Results have been obtained for the shapes of the 
bubbles of different sizes, at different depths of the bath. The 
transient concentration profiles both inside and in the mass 
boundary layer have also been obtained. The variations of drag 
coefficient with Reynolds number and EBtvfis number were also 
calculated with bubble rise. The results are presented and 
discussed below. 


4. 1. Validation of the Model for shape cal culations 

In fig. 4.1, the shape of a falling rain drop has been 
presented. Experimental data for the shape of bubbles in liquid 
metal bath were not available in literature. This because of the 
formidable difficulties associated with taking photographs in 
opaque and hot liquid metal medium. Therefore to check the 
accuracy and validity of the computational scheme used for the 
prediction of shape, the shape of a falling drop in air was 
predicted.. Figure 4.1b is the numerically predicted shape which 
matches to a very good accuracy with the experimental photograph 
of the rain drop of the same size, provided by Choji and Morogon 
C53Mac3. From the comparison, it is evident that the 
computational scheme we have used for the shape calculation is 


correct. 



FIG. 4.1a. ACTUAL PHOTOGRAPH OF RAINDROP 
(DIA. 6 mm) [PROVIDED BY CHOJI 
MAGONO] 



FIG. 4.1b. NUMERICALLY PREDICTED SHAPE OF 
RAINDROP OF SAME SIZE. 



Depth in meter 






FIG. 4.2. SHAPE OF RISING BUBBLE AT DIFFERENT DEPTHS. 





4.2 Results and discussions . 

Figure 4.2a to c , present the shapes of rising bubbles 

in a liquid metal bath. From the figures it is clear that as the 

bubble size increases, the deformation becomes more. Since the 

static pressure decreases considerably with bubble rise, the 

volume of the bubble increases. Consequently the deformation 

becomes more severe. The bubble deforms from spherical to obloid 

ellipsoidal shape. When the bubble reaches the free surface the 

volume drastically changes. For this reason it has not been 

possible to predict the shape of the bubble very much near to the 

free surface. In general the bubble has flat bottom and a curved 

top surface. The flat bottom surface is believed to be a 

consequence of the dynamic pressure ( in liquid) being very low in 

the wake behind the bubble. This explanation is corroborated by 

the streamline patterns discussed later. With respect to the 

initial value of equivalent bubble diameter (d ), it is seen that 

eq 

larger the value of d^, larger is the departure form spherical 
shape. This is not surprising as larger bubbles experience 
greater buoyancy and inertial effects. The increase in buoyancy is 
only due to size increase while that of inertia is both to the 
size and velocity change. 

Figure 4.3 shows the variation of the drag 
coefficient with Reynolds number, for different EStvtts numbers. As 
the Reynolds number increases , the drag coefficient decreases. 
With respect to Ettttt vs number on the other hand Cp increases. The 
effect of Reynolds number on Cp is well known. For low Reynolds 
number Cp varies as 1/Re and the present study indicates that this 
variation is approximately true for bubble also in the range (K Re 
^ 100. The drag coefficient decreases with Re because of the 
weakening of viscous forces in comparison with inertia, at higher 
Re. The drag coefficient increases with ES because of the larger 
contribution from the form component of the drag force . The 
deformation becomes more at higher ES and hence the form drag 



A = 0.0091 
B = 0.0094 
C = 0.0099 
D = 0.019 
E = 0.038 
F = 0.1 
A':- 0.009 
B'=-0.013 
C's- 0.052 
D =-0.081 
E- = - 0.092 
F =-0.132 


FIG. 4.4a. STREAMLINES AT 

Re outside s 100, Re inside i 4.0, 
Edtovs number = 3.12 




FIG. 4.4d. STREAMLINES AT 

Re outside =100, Re inside = 4.0, 
Eotovs number = 1.176 



Cu 2 S 



FIG. 4.5 DISTRIBUTION OF Cu 2 S IN MASS TRANSFER 
BOUNDARY LAYER. 



component increases. 

Figure 4.4 a to d presents the stream lines both 
inside and outside the bubble for different Reynolds numbers. It 
is seen that the flow field and the shape are very sensitive to 
the Reynolds number. The inside flow pattern changes with shape to 
a large extent. However the vortex center is more or less located 
in the vicinity of 0.7 (non-dimensional radius), as for a 
potential flow vortex inside the bubble. At the external Reynolds 
numbers of 20 and 50 , the circulation pattern inside the bubble 
exhibits fore and aft symmetry. At Re Q = 100 , the vortex moves 
closer to the front stagnation point, possibly due to maximum 
surface shear sources in the front region, at higher Re. This 
asymmetry may be due to the separation of external flow. As the 
external Reynolds number crosses about 50, the flow separates 
behind the bubble. At high Ett , however, the critical Reynolds 
number for flow separation decreases due to larger bubble 
deformation. For solid spheres, it is known that flow separation 
occurs at a Reynolds number of around 20. For bubbles, this limit 
is higher because they offer no viscous tangential stress. Beyond 
Reynolds number 100 and ES number about 3 , the shape of the 
bubble becomes instable. For this reason it was not possible to 
predict the spherical cap bubbles which are generally encountered 
in metallic systems. 

Figure 4.5 shows the variation of concentration of Cu^S 
with angular positions within the mass transfer boundary layer. It 
is observed that there is a sharp boundary layer in the frontal 
regions of the bubble, whereas the boundary layer is very thick 
near the rear stagnation point. These findings are in conformity 
with those of Chao C69Cha3. 

Figure 4.5a shows the variation of oxygen concentration 
along & = n/Z at 0.04 nondimensional time for different gas phase 
®ass peclet numbers. It is clear that the maximum concentration 
happens to be near the vortex center, this being a stagnation 
zone. Again with the increment of Re.^ and Pe i this zone becomes 
narrower and sharper. The profile indicates the existence of a 


DIMENSIONLESS CONCENTRATION OF 



RADIAL DISTANCE 

FIG. 4.5a. 0 2 CONCENTRATION IN GAS PHASE AT 6 - 7T/2 

ALONG RADIAL DISTANCE. 




FIG. 4.5b. CONCENTRATION PROFILE AT D 
LOCATION AT FIXED TIME . 




FIG. 4.5c. VARIATION OF CONCENTRATION PROFILE 
AT 7T/2 WITH TIME. 




FIG. 4.5ci 


VARIATION OF CONCENTRATION AT ONE POINT 
WITH TIME. 



FRACTION OF 0 2 REMAINING 



FIG. 4.5e. FRACTION OF 0 2 LEFT IN GAS PHASE AT 
DIFFERENT TIME. 



teep mass transfer boundary layer near the bubble surface, an 
xygen deficient internal wake near the axis and an oxygen rich 
ore. These features have been theoretically predicted and 
iscussed in detail by Brignell C75Bri3. Figure 4.5b presents the 
ariation of oxygen at different B along the radial distance. It 
s observed that as B nears the front stagnation point the 
egion is depleted of oxygen . This is because of the strong 
irculation that is established inside the bubble. Thus the air 
hat comes from the frontal region to the rear region is devoid of 
xygen due to reaction at the interface . Figure 4.5c shows the 
epletion of oxygen with time, at & = n/Z along the radial 

istance. Here as time progresses the radial distance at & = n/Z 

ncreases due to deformation. However the concentration at the 
ortex center remain at maximum for all the time steps. Figure 
.5d presents the variation of 0^ at a point at different times, 
t is noticed that initially the concentration falls down rapidly 
ut as time progresses the drop becomes slow. This is because 
nitially depletion rate of oxygen is very high since the amount 
f oxygen available is high. But towards the end, the reaction 
lows down due to lack of oxygen inside the bubble. For lower 
nside mass peclet number the depletion of oxygen is faster 
ecause of higher rate mass transfer. The above mentioned ideas 
re reinforced in the figure 4.5e which shows the variation of 
vailable oxygen fraction within the bubble with time . 

The transport , flow and shape change results predi cted 
In the present work, thus are in conformity with the observations 
earlier investigators . Elaborate comparison however, have not 
been possible due to the paucity of detailed results in the 
literature . The present study thus off ers a powerful solution 
technique for solving the complex flow and transport processes in 
and around non— spher i cal bubbles. 



CHAPTER-V 


CONCLUSIONS AND SUGGESTION FOR FUTURE WORK 

From the results obtained in the present study, the 
Lng conclusions can be drawn. 

- The present model offers a powerful solution technique 
Lving complex flow and transport processes in and around 
lerical bubbles for flows in low Reynolds number range (O^ 
). 

- The shape of the bubble , and flow, in and around are 

Lely connected. The shape is sensitive to Reynolds number 
irge extent in the range 0<Re£ 50. Beyond this the 

ition is not large . For Re of the order of 200 the scheme 
L o give a converged solution of the flow. These two aspects 
so dependent on E&tv&s number . 

—The shape changes considerably near the free surface of 
Lh. This is because of the sudden increment in the volume of 
able . 

— The basic nature of the transport is not signif i cantly 
?nt from that of the spherical cases studied by earl ier 
i except in order of magnitude 
Lions: 

The present study assumes that the flow is laminar . But 
ial situations this may not be true . To simulate the 
:al situation it may be necessary to take the turbulent 
into account. 

High Reynolds number should be tried. For this modified 
5 e . g . upwind FEM schemes can be tried to obtain converged 
an to the flow. 

Here a single bubble behavior has been studied. Multi- 
system can be studied to simulate practical situations . 

A better interpolation technique can be developed which 
restricted to small deformations . 
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APPENDIX -I 


DEAN CURVATURE OF A SURFACE OBTAINED BY REVOLVING A CURVE 
AN AXIS : 

Consder the surface obtained by revolving a curve in 
X Z ) plane about Z-axis. Denoting by & and 4> the polar and 
thal angles respectively, the equation for this surface in 
ical co-ordinates is given by r ■ R C6 ), where R i & ) is a 
function of 0. The mean curvature of a surface at a point P 
is defined as 

KCP5-IK C n.5 + K C 1 Cl!) 

mi m2 

K ( n . > and K <n„) are the normal curvatures along two 
n 1 n E ^ 

ly perpendicular directions n^ and n^ tangential to the 
:e at P . The normal curvature at a point on the surface along 

A 

rection n , tangential to it at that point is given by 

K Cm3 -C^fb.m C23 

m n 

N is the unit normal to the surface at the point and — 

3 

> R = E n .V . T5 is the rate at which the unit normal to the 
i=1 1 1 

e turns as it is moved along the direction n . It can be 

a 

that the normal curvature K^Cn) along a direction is the 
.ire of the curve on the surface obtained by its intersection 
le plane containing TQ and n . It can also be shown that 
jh only two mutual ly perpendicular di rections n ^and n ^ on 
igent plane are used to evaluate the mean curvature, the 
ilue will be obtained when the mean is cal culated over all 
ons tangential to the plane at the point. In this sense it 
y the mean curvature for it is the mean value of the 
re of al 1 normal curves C i. e . curves obtained by the 
ction of the surface with the plane containing the normal 
surface ) on the surface and passing through the given 
We now proceed to calculate the mean curvature at a point < 

) on the surface, defined in eqn-1 . At such a point the unit 
R is given by 



? / ( > 
I 9 / I 


1 = 


( r, S, )= 0 is the equation of surface . In the present 
have / ( r, T 4> > = r- R ( © ) , so that 

= * wa f n p A /p 

R = ( r - R/R d >/( 1 +R C /R C ) ,/c 

* * A p *p W /p 

= (R r - R ©)/ <R + R ) .7 <3) 

9 = r &/&r + a/r &/&& + <£/r sin a &/&4> in the spherical 
tes, r, a, arethe unit vectors in the di rections of 
ng r , d and 4> respectively. In terms of the unit vector i » 
along the three cartesian axes x r y and z, we have 


r = sin £ cos <jS> i + sin a sin j + cos a k 

^ «*W «IV alV 

a = cos 9 cos i + cos 0 sin ^ j - sin a k 

4 , = - sin i + cos 4 > k 


(4) 


mutually perpendi cular unit vectors n ^ and n^ in 
1 be chosen as: 

n z = 4> and n 1 = n z >M = CR a + Rr>/<R 2 +R' 2 ) 1/2 

(5) 

lation of the normal curvatures along n^ and n^ we need 
ttors n„j. 9 and n^.^, which are easily obtained as 

,.5 = i/(r z + r , 2 , 1/z [ r ' *'*' + •'» B ] <6> 

9 = 1/ R sin a S/S<k 

be noted that in the above expressions for ^.9 and 
e have put r = R ( a ) , corresponding to the interface . 
(3) and (6) 


= 1/CR 2 +R 

2 *2 
1/(R +R > 


*2 1/2 * ' * p • o -d/o 

) a/*a (R r - R a)/ CR e + R ) c 

[ (Rr+R d r/& & - R "& - Rd &/&&)/ (R 2 +R* 

" » * • « _ 2 ’2 3/2 

— (R r - R a>(RR +R R l/CR^+R 





i eqn <4) we have 


S r/&& = © and & &/&& = -r so that 


9 )R = 1 / { R+R 


* 2 )|^2 R* 


r + ( R-R )© -CR r-R6)(RR+RR ) 

2 *P 

CR c + R c ) 


= 1/(R 2 + R* 2 ) 2 |\r 2 +r' 2 ) (2Rr + <R-r")© - ( RR*+rV’><. 


R r - R & ) 


p p $ * *i * p *9 O «» *1 

1/CR+R c ) (R K + 2R -R R R ) r +( R +2R R -R C R >© 

— = P *P 1/PF * • * 1 

(n 2 -9 ) N = 1/Rsin©. 1/CR+R c ) R S r/&4> - R a &/&$ ..(7) 


ing S r/&4> = sin© 4> and & &/&<£ — cos© we obtain 


9 )R = 1/Rsin© .1/CR 2 +r’ 2 } 1/2 f( Rsin© - fi'cos©) 


illy 

* — - — p • jj c /p r p ®p « a * «i 4 * 

Cn „ ) = n„- [(n,. 9 >W □ = 1/<R C +R ) . R+R +2R -RR R +R + 

» 1 11 L p .p o .. 

+2R R -RR 


p *P *5 /Pf P *P P P «p • A *p “ o I 

= 1/CR+R ) CR +R )+RR +R -RR R — R R 


p •? s/?r p *p p *? p *p *■ p * p 
= 1/CR+R )! (R+R )+R CR+R }-RR (R^+R ) 


2 *2 5/Pf P '? ' 

1/CR+R )! R+2R -RR 


(n_) = n_. C ( n_. v m 2 = 

i 2 2 2 


Rsin© - R cos© 

2 *2 1/2 
R sin© C R+R > 


t stands K (n„> can either be positive or negative depending 

^ ^ f m 

i whether Rsin© > R cos© or R sin© < R cos© . However if ^ is 

outward normal to the surface which is true of H in eqn (3> , 

i K (n„) is always +ve . Hence we take the modulus of the RHS 
n 2 

:onsidered in the expression for K (n„) . Thus , 


Tig- C ( Tig. 9 >N □ 


j Rsin© - R cos3 


a w /p 

R sin© ( R c +R ) 


( 9 ) 


in© is always +ve for 0 < © < n ) 
obtain 

2 »2 
+2R -RR 

r 2 + r' 2 > 3/2 

an curvature at a point on the surface obtained 
the curve r = R (© ) on the X— Z plane about Z— axis. 


| Rcos© - Rsin© J 

2 *2 -i /2 

R sin© ( R +R 


( 10 ) 


by 



APPENDIX- 2 


>UPLING OF ENERGY EQUATION WITH MOMENTUM EQUATION : 

Using Boussinesq approximation one can proceed with the 
sumption that the density varies in the buoyancy term of the 
imentum equation and every where else there is a average 
nsity p. 

en pressure can be expressed as 
^ ^static + ^dynamic 


ere static pressure is 

P is the average density of the fluid. 


Static = P 0 + ^ 9 z 


p„ is atmospheric pressure 


w 


-Vp 


{ p q — p q )k - V p 


dynami c 


= ( P - p! g k - 7 p 


dynamic 


ange in density A p can be expressed in terms of coefficient of 
ermal expansion as fallows : 

ft = - 


dV 


V dT 

i . e. , change in volume per unit* volume per unit temperature 
Ff erence ) 

Where V = — - 


<i ft = p ( 


1 dp 

2 dT 
P 


1 dp 


dT 


or { p — p) = p f 3 ( T - T ) 

>re T is the reference temperature 


sn, 


-Vp = p q ft ( T 


T . )k 
ref 


V p 


dynam i c 


t coupling momentum equation with the energy equation we have 


<7 . V = 0 ( 1 > 

— + P« V .V V ) = P g fi ( T- T ref )k - V P dynafflic + 

+ 7 . ( p (7 V V T ) ) (2) 

+< V . V T > = <x ( ? 2 T ) + S -(3) 



APPENDIX-III 


CONDITION AT INFINITY OF OSEEN TYPE 


One way of reducing the flow domain without losing accuracy 

to consider an Oseen type of flow at large r- As suggested by 

gers C85Chr3 the flow at infinity can be thought of being 

ated from a point force applied to the origin of the coordinate 

tern. The magnitude of this force is equal to the drag force 

rted by the liquid on the bubble . 

At suf f i ciently large distances from the bubble, the flow 

scity differs slightly from the uniform profile U . The 

00 

jarised Navier— Stokes equation regardless of the magnitude of 
Reynolds number is 


C cos & — — 
Sr 


sin & S w _ 2. «. 

—r~ >f " " D < 


Appropriate solution with the required symmetry is given by 


= -HC„ (cos a + 1 ) El— exp£-rU ( 1-cos a>>! <H) 

2 2 ce> 

LbWa.nl 


This solution gives a flow far from an arbitrary three 
nsional body with no lifting force. In the present case this 
ies because the lifting force is absent owing to the symmetry 
he flow. The constant depends only on certain details of the 
in the vicinity of the body, namely the drag force. The later 
upposed to be equal to the buoancy force . 


F a = ip l~ p g > ^ 

b pg is the gas density which will be neglected as 
ared to the fluid density p^« Also V denotes the volume of the 
Le. 

Now assuming that a sphere of radius fa is situated at the 

in of the coordinate system and b is selected so that the drag 

i exerted from the liquid on the bubble is equal to the drag 

i exerted on this imaginary rigid sphere , then 

F = 6 nfjhU 
a oo 



or, b 


(4) 


= __2V 

isnv U 

Ob 

It remains only to match the constant Cg with this 
solution. For small r the stream function from (2) is represented 
as 

¥>2 = -H Cg ( 1+cos & ) rU^ /(2i>) < 1-cos & )= -Cg rU^ fv sin 2 & <5) 

This solution has to be matched to Stoke's solution f or a 
sphere of radius b, which is given by 

2 

¥> = -3/4 rbU sin & 

00 

From the above expression it is immediately seen that 

Cg = 3/4 

and thus 

2 2 

V* = 1/2 U r sin & — 3/2 ub( 1+cos & ) C1-exp£-rU A* (1— cos 

From (4) = 1/2 U r 2 sin 2 & - - (1+cos 0) El-expC— rU / v 

00 (1-cos &)}1 (6) 

The velocities can be determined from the definition of 
stream function at the cut off region and assigned as the boundary 
conditions. 



APPENDIX-4- 


frontal METHODS 

When the number of unknowns in th. 

large, the method adopted for solvino ,e Pr0Me "' is 

significant bearing on the computer ^ 

execution time. It would be uneconomicai to sol ^“ lri "" ent and 

hy simple matrix inversion technigues. The Frontal 

on direct gaussian elimination method for i""' 1 '' 0 ' 1 “ b *” d 

matrices where the leading diagonal is always used ' 

unsymmetric matrices, encountered in a „i„ 35 PiVOt - For 

problems , the most suitable p ivot • ra " 9e ° f smearing 

leading diagonal; and frontal s 0 i uti n0t , neCI!S5ariI y on the 

diagonal pivoting CSITayU. This h Q " r ° Utlne exist s for off 

y ,nxs - however, tends tn t. 

consuming . The method incorporated here use "" ““ 

pivoting , but includes many feature, of the ’ °" Iy di * alm * 1 
solving unsymmetric matrices. ' r ° n al tech nique for 


In general terms , the overall c n i * ■ 

(1) Fo nation of Elemental Matrices. U 100 tSChnique consists of 
Assembling into a global matrix. 

<3> Intr0ducti °n of boundary conditions. 

C4) Reduction of global matrix usinn 

technique. 9 Gaussian elimination 

<5) A back substitution process. 


1BSIC P HILOSOPHY OF THE FRQNTAI METHOD ; 

The primary objective of the fm^t ^ 

elimination of variables as soon a * method 15 the 

soon as possible 

Introduction into the global matrix via the .on • their 

When the contributions from all the elements to 

point have been assembled, the corresnona “ l * r " od *l 

corresponding variables 

With that node can be eliminated. The complete matrix is 

“ ■ ln “ " dd ‘ ad equation can be .Hi” t.™ 

e core and stored in the disc. The equations held i„ core "" 

the corresponding node, and variables are collectively J led 

fr °"t 3 " d the »—■ ’ df unknown variables in th. front “ l. 1 

M fr “ nt “ idth - Thb a continuaily . 



O Inactive Nodes 
• Active Nodes 
© Deactivated Nodes 

^3 Front 


Vs '/A Next 


Element for Assem 



Assembled 


Element 


bly 


Fig. A4. Definition: of Front and Element Numb 
for Minimum Front Width. 



after the contribution to a particular node have been fully 
summed, the corresponding equation reduction based on a diagonal 

pivot can be executed. 

A pre assigned global matrix core area is first filled from 
contributing elements and the largest diagonal entry in the pre 
assigned core is found and used as pivot in a direct gaussian 
el iminat ion process . As the maximum predetermined number of 
equations are eliminated, the corresponding reduced equations are 
written on to the disc and more elements and corresponding 
equations are taken into the core. Fig-A4 shows the basic idea of 
the frontal method. The equations , nodes , and variables currently 
in the core are termed “active " ? those assigned to the disc are 
"deactivated" ; and those yet to appear in core are "inactive" 
E73Tay3 . 


/ r\ 


APPENDIX - 5 


TRANSFINITE INTERPOLATION: 

In two dimensions a linear Lagrange 
interpolation function can be written in terms of 
interpolations in each curvilinear direction as: 

2 

£< K » n > = T ^ (-y— ) r ( ? , t>) (1) 

n^1 



2 

r-l £ , = y v ( r ( ? , y> > 

— 4dA 5ft J ““ 111 


C2) 



This interpolation is called transf inite since it matches the 
function on the entire boundary defined by ? = 0 and if = I in the 
first equation or by rj = 0 and r> = J in the second. 

The tensor product form is 


2 2 


n > 

= Y Y 

n— 1 4=1 

V 

>nm 

*n ’ 

V 





0 


> r 


( 3 ) 



The sum of eqn(l) and eqn(2) is 


§( Jf, Y>) 



) r(£ 


n 


rf}) 


2 





) r ( 




nj 

m 


( 4 ) 


Equation (4), when evaluated on £ =0 boundary gives 

2 

S ( 0. y> > = r ( 0, y>) + y y C-4~> r( Q, n ) (5) 

— — l, ns J — in 

m-1 

This does not match the function on the £ = 0 boundary because of 
the second term in the right, which is an interpolation between 
the ends of this boundary 


r( 0, J ) 




r( 0, r>) > < T v (-r >r < 0, » ) 

— Z* , m J — m 

m-1 


'?0,1) 


Similar effects occur on the other boundaries and the discrepancy 

on £ = I boundary is 2 

” " ) 


7 V (— -t — ) r< I, » ) 

4=1 m J 


The discrepancies on both these boundaries can be removed by 
subtracting from S( Z , rj ) , a function formed by interpolating the 
discrepancies between the two boundaries. 


Ri Z , r,) 


= y t %— > y v < — a > r( ? » y) ) 

n 4l n 1 L m 4 1 m J - n ra 


r 2 

r v 


( 6 ) 


This is simply the tensor product form given by eqn (3) 
which matches the function at the four corners. 


The function S-R then matches the function on all four sides 

mm m m 

of the boundary, so that we have the t ransf inite interpolation 
form. 


Z, r>) = 


!,V 


n— 1 


t->£ < V 


Ji L V 


n> + l/n'-T 

m-1 


ri 


■> 


> ¥* t— r-> r C Z . Y> ) 
Hi J — n m 


(7) 


which matches the function on entire boundary. 
The tensor product form 


r(f , Y) > = 


n 


2 2 
=4 m— 1 


n 


¥»< 

m 


Y) 


> r ( Z 


J ' - 


n 


n 


m 


( 8 ) 


matches the function at the four corners of the boundary. This 
generalises to the interpolation from a set of N+M intersecting 
curves for which the univariate interpolation is given by 


r 


N 

,Y>>=E * 

1 n 



(?„ ,Y) ) 


n 


( 9 ) 


M 


2 



( 10 ) 


- (1: ' n>= E' B m E -"no 


where now the blending functions, <p and are any functions 

Tl ill 

which satisfy the cardinality conditions 

4> (“1 )=■$- n=1 ,2, N 1=1,2, N 

Mi I nl 

*m (~ )=£ _ m=1,2, M 1=1,2, M (11) 

J mi 

The general form of the transf inite interpolation then is 
N M 

E «T ■«>=£ * n <-f-> - «„ .« > + E (-?-) E <JT.n„ > 


or. 


N M 

~ E E V_ (- 4 -» r (* ,Y) ) 

n=1 m=1 n 1 m J - n m 

M N 

i <f ■« >=E *•„ (-3-) £ <« .*>„) + E *„ <- r > 
[ e «„■»>- t *•„ <- 3 -> c )] 

L m=1 J 


( 12 ) 


The first term is the result at each point in the field of 
the unidirectional interpolation in the r> direction and the 
bracket is the result of the difference between the specified 
values on the ?=£ lines and result of the unidirectional 

O 

interpolation on those lines. The two directional transfinite 
interpolations can thus be implemented in one di rect ion , say r?, 
over the entire field. Calling the result 


v<’»>=£ *v <-3-> - « -v 

- 1 


(13) 


and interpolating the discrepancy on the lines over 

entire field in the other direction, with the definition 


N 

F a (*,*)=£*_ 
- 1 n 
we get 


(- 4 -) E> 

I — 


<? n ^ )_F i ( * n 


(14) 


the 


r (e,^)^ (?,Y))4F 2 (?,r>) (15) 

This transfinite interpolation form is the best algebraically 
approximated C85ThoU. 
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